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Scattering hinders the passage of light through random media and conse- 
quently limits the usefulness of optical techniques for sensing and imaging. 
Thus, methods for increasing the transmission of light through such random 
media are of interest. 

Against this backdrop, Dorokhov [1], Pendry [2,3] and others [4,5] theo- 
retically predicted the existence of a few highly transmitting eigen-wavefronts 
with transmission coefficients close to one in strongly backscattering random 
media. The breakthrough experiments of Vellekoop and Mosk [6,7] confirmed 
the existence of these highly transmitting eigen-wavefronts and demonstrated 
that they could be discovered by using information from the far side of the 
scattering medium. 

Here, we numerically analyze this phenomenon in 2-D with fully spectrally 
accurate simulators and provide rigorous numerical evidence confirming the 
existence of these highly transmitting eigen-wavefronts in random media com- 
posed of hundreds of thousands of non-absorbing scatterers. 

We then develop physically realizable algorithms for increasing the trans- 
mission through random media using backscatter analysis. We show via numer- 
ical simulations that the algorithms converge rapidly, yielding a near-optimum 
wavefront in just a few iterations. We also develop an algorithm that combines 
the knowledge of these highly transmitting eigen-wavefronts obtained from 
backscatter analysis with intensity measurements at a point to produce a fo- 
cus using significantly fewer measurements. 
© 2013 Optical Society of America 
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1. Introduction 

Media such as glass and air are transparent because light propagates through them without 
being scattered or absorbed. In contrast, materials such as turbid water, white paint, and egg 
shells are opaque because the randomly arranged particles cause light to scatter in random 
directions, thereby hindering its passage. As the thickness of a slab of highly scattering 
random medium increases, this effect becomes more pronounced, and less and less of a 
normally incident light is transmitted through [8]. 

In this context, the theoretical work of Dorokhov [1], Pendry [2,3], and others [4,5] pro- 
vides unexpected insight into how, and the extent to which, the limitations imposed by 
random scattering may be overcome. Specifically, these authors predict that in highly scat- 
tering random media composed of non-absorbing scatterers, the eigen-wavefronts associated 
with the right singular vectors of the S21 or transmission matrix will have transmission co- 
efficients whose distribution has a bimodal shape as in Fig. 2. Consequently, while many 
eigen-wavefronts have a small transmission coefficient, a small number of eigen-wavefronts 
exist that have a transmission coefficient close to one, i.e., they propagate with almost no 
scattering loss. 

The breakthrough experiments of Vellekoop and Mosk [6,7] provide evidence of the exis- 
tence of these highly transmitting eigen-wavefronts in random media. Vellekoop and Mosk 
showed [6] that intensity measurements on the transmission side of a scattering medium 
could be used to construct a wavefront that produced about 1000 x intensity enhancement 
at a target point over that due to a normally incident wavefront. Their work set off a flurry 
of research on experimental methods for measuring the transmission matrix and compar- 
ing the transmission coefficient distribution with the theoretical prediction [9-13], faster 
experimental methods for focusing [14-18], and numerical work on the properties of the 
eigen-wavefronts [19]. 

Our work is inspired by these three lines of inquiry. We develop iterative, physically realiz- 
able algorithms for transmission maximization that utilize backscatter analysis to produce a 
highly transmitting wavefront in just a few iterations. These algorithms build on the initial 
algorithm presented in [20]. Furthermore, we develop an iterative, physically realizable algo- 
rithm for focusing that utilizes intensity measurements at the desired point and backscatter 
analysis to produce a wavefront with significantly fewer measurements than other approaches. 
A crucial feature of the algorithms we have developed is that it allows the number of modes 
being controlled via a spatial light modular (SLM) in experiments to be increased without 
increasing the number of measurements that have to be made. 

Finally, we numerically analyze the phenomenon using a spectrally accurate simulator for 
2D scattering systems and provide the first numerically rigorous confirmation of the shape 
of the transmission coefficient distribution and the existence [7] of an eigen-wavefront with 
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transmission coefficient approaching one for random media with a large number of scatterers. 

The paper is organized as follows. We describe our setup in Section 2. We discuss the 
problem of transmission maximization and focusing in Section 3. To assist in the develop- 
ment of physically realizable algorithms for these applications, we identify physically real- 
izable operations in Section 4, and describe iterative, implementable algorithms for finding 
transmission-maximizing and focusing inputs in Sections 5 and 6, respectively. We highlight 
the existence of the eigen-wavefronts with transmission coefficients approach one, the al- 
gorithms' performance and rapid convergence via numerical simulations in Section 7, and 
summarize our findings in Section 8. 

2. Setup 

We study scattering from a two-dimensional (2D) random slab of thickness L and periodicity 
D; the slab's unit cell occupies the space < x < D and < y < L (Fig. 1). The slab contains 
N c infinite and ^-invariant circular cylinders of radius r that are placed randomly within the 
cell and assumed either perfect electrically conducting (PEC) or dielectric with refractive 
index n^; care is taken to ensure the cylinders do not overlap. Fields are TM^ polarized: 
electric fields in the y < (i — 1) and y > L (i = 2) halfspaces are denoted e^p) = ej(p)z. 
The field (complex) amplitude ej(p) can be decomposed in terms of +y and — y propagating 
waves as ej(p) = ef(p) + e~(p), where 



/v 



ef(p) = £ h 



n u in c 



-N 



xx + yy — (x,y), k n — k ntX x i k nt yy — {kn,xi i^n,t/) ; k n 



\kn\\2/k nj y is a power- 



In the above expression, p 

27rn/D } k nty = 2ir^(l/X) 2 — (n/D) 2 , X is the wavelength, and h n 
normalizing coefficient. We assume N = \_D/X\, i.e., we only model propagating waves and 
denote M = 2N + 1. The modal coefficients af n , i = 1, 2; n = —N, . . . , N are related by the 
scattering matrix 

•Sn 5*12 
S21 S22 



(2) 



where af = af_ N . . . af . . . af N . In what follows, we assume that the slab is only 
excited from the y < halfspace; hence, = 0. For a given incident field amplitude ef(p), 
we define transmission and reflection coefficients as 



15*21 ' <±i\\2 



(3) 



3 



and 



l^n ' ^i Hi 



r(o+) := " 1 + - 1 2 " , (4) 

NhI 1 1 2 

respectively. We denote the transmission coefficient of a normally incident wavefront by 

i T 



^normal 7~\ 



f 



); here T denotes transposition. 



3. Problem formulation 

3. A. Transmission maximization 

The problem of designing an incident wavefront a opt that maximizes the transmitted power 
can be stated as 

a opt = arg max r(af) = argmax ^ ^ + ~2 ~ = argmax \\S21 • afWl (5) 



_1 1 



H«l 112 ||a+[| 2 =l 



where || \\2= 1 represents the incident power constraint. 

Let 5*2i = Yl!iLi a iM.i ' V4 denote the singular value decomposition (SVD) of S 2 i; o~i is the 
singular value associated with the left and right singular vectors and v i: respectively. By 
convention, the singular values are arranged so that o~i > . . . > &m an d denotes complex 
conjugate transpose. A well-known result in matrix analysis [21] states that 

a opt = JLv ( 6 ) 

When the optimal wavefront a,, t is excited, the optimal transmitted power is r opt := r(a t ) = 
o~\. When the wavefront associated with the i-th right singular vector is transmitted, the 
transmitted power is t{v?) = erf, which we refer to as the transmission coefficient of the i-th 
eigen- wavefront of ^21 . Analogously, we refer to r(uj) as the reflection coefficient of the i-th 
eigen- wavefront of S 2 i- 

The theoretical distribution [1-5] of the transmission coefficients for lossless random media 
has density given by 

1 M /I 
/(t)= hm J2 5 (r-r(v i )) = — —==, for 4exp(-L/2/) < r < 1. (7) 

i=l 

In Eq. (7), I is the mean-free path through the medium. Fig. 2 shows the theoretical density 
when L/l = 3. From, Eq. (7) we expect r opt = 1. 

From (6) it follows that the optimal wavefront can be constructed by measuring the S21 
matrix and computing its SVD. Techniques for measuring the S21 matrix have been developed 
in recent works by Popoff et al. [9,10] and others [11,12]. Kim et al. experimentally measured 
the S21 matrix and demonstrated improved transmission by using the optimal wavefront in 
Eq. (6) [13]. 
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In the lossless setting, the scattering matrix S in Eq. (2) will be unitary, i.e., S H ■ S = I, 
where I is the identity matrix. Consequently, we have that Sn ■ Sn + S$[ ■ S 2 i = I, and the 
optimization problem in Eq. (5) can be reformulated as 

Oopt = argmax (af) H ■ S 21 ■ S 2 i • at = argmin \\Sn ■ af\\ 2 = argminr(af ). (8) 

ll^l|2=l ' ' " " / ||a||| 2 =l a+ 

=(a+)«-(/-Sf 1 .S 11 ).a+ 

In other words, in a lossless medium the backscatter-minimizing wavefront also maximizes 
transmission. Let Sn = 'V.f denote the SVD of Sn; a.j is the singular value 

associated with the left and right singular vectors and v i , respectively. Then from [21] it 
follows that 

a op t = 2m- (9) 

When this optimal wavefront is excited and the medium is lossless, r opt = 1 — r(a opt ) = 
1 — a M = o\. When the wavefront associated with the i-th right singular vector v_ { is excited, 
the transmitted power is given by = 1 — r(^) = 1 — af, which we refer to as the 

transmission coefficient of the z'-th eigen-wavefront of Sn. Analogously, we refer to Tfy) as 
the reflection coefficient of the z'-th eigen-wavefront of Sn- 

A technique for increasing transmission via backscatter analysis would require measure- 
ment of the Sn matrix and the computation of a opt as in Eq. (9). Our objective is to develop 
fast, physically realizable, iterative algorithms that converge to a opt by utilizing significantly 
fewer backscatter field measurements than the O(M) measurements it would take to first 
estimate Sn and then compute its SVD to determine v_ M . 

3.B. Focusing 

From Eq. (1) and using the fact that that a 2 = S21 ■ (since — 0), the field at point p Q 
is 

4(P ) = \h_ N e- jkt ^o ■ ■ ■ h N e- jk+ ^o] -S 21 • a+. (10) 
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The problem of designing an incident wavefront that maximizes the the intensity (or ampli- 
tude squared) of the field at p Q is equivalent to the problem 



Ofoc = argmax— — ^— = argmax ||/ (pj ■ S 2 i -a^|||, (11) 

at WSLl lb 




whose solution is 



a ^ ~ S?1 ' L ^ ] (12) 

- foc ~ mp )\\2~ m-i(p r 
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Thus the optimal wavefront equals the vector c(p ) with normalization to satisfy the 
power constraint. It can be shown that this wavefront may be obtained by time-reversing 
the wavefront received by placing a source at p Q [22]. This fact was exploited in recent work 
by Cui and collaborators [23,24]. 

In Vellekoop and Mosk's breakthrough work [6,7,25], a coordinate descent method was 
employed for constructing the optimal wavefront. The coordinate descent approach finds the 
amplitude and phase of a single mode that maximize the intensity at p Q while keeping the 
amplitudes and phases of the other modes fixed and then repeating this procedure for the 
remaining modes, one mode at a time. In Vellekoop and Mosk's experiments [6,7,25], they 
kept the amplitude constant for all the modes and considered phase-only modifications of 
the incident wavefront. While this reduces the complexity of the algorithm, this approach 
still requires O(M) intensity measurements at p Q to construct the optimal wavefront. When 
M is large, the time for convergence will also be large. 

This has motivated recent work [16-18] for faster determination of the optimal wavefront. 
Cui [16, 17] considers an approach using multiple frequencies to find the optimal phases 
of modes simultaneously, while Stockbridge et al. [18] have proposed a coordinate descent 
approach using 2D Walsh functions as a basis set. These methods have accelerated the 
experimental convergence, but the reported results are still for small M (between 441 and 
1024). 

Expressing the optimal wavefront in terms of the singular vectors of S 2 i yields the expres- 
sion 

M M 

flfoc oc S£ • /(p Q ) = a i fef • /(Pq)) m« = Yl WOk- ( 13 ) 




i=l 



Recall that erf = t(i^); thus an important insight from Eq. (7) and Fig. 2 is that most 
of the singular values in Eq. (13) are close to zero. However, there typically are K <C M 
singular values close to one. It is the superposition of these K eigen-wavefronts of S%i having 
transmission coefficients close to one whose constructive interference yields the maximal 
transmission that contributes to maximal intensity. 

In the lossless setting, when the scattering matrix S is unitary, we have that r{v_i) — 
1 —T(v_M-i+i)- Hence, the K eigen-wavefronts of S21 that have transmission coefficients close 
to one correspond precisely to the K eigen-wavefronts associated with Su that have reflection 
coefficients close to zero. By using O(K) backscatter field measurements to measure the K 
eigen-wavefronts of Su with small reflection coefficients and 0(K) intensity measurements 
at p , we might expect to approximate a foc in Eq. (13) and yield a near-optimal focus using 
just O(K) measurements (we expect K <C M). 

Our objective is to develop a fast, physically realizable, iterative algorithm that utilizes 
backscatter field measurements and intensity measurements at p to construct a near-optimal 
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focusing wavefront using significantly fewer measurements than are required by coordinate 
descent methods that only employ intensity measurements at p Q . 

4. Recognizing physically realizable matrix-vector operations 

The iterative algorithms we will develop in Sections 5 and 6 build on the vast literature 
of iterative methods in numerical linear algebra [26,27]. The algorithms are based on three 
matrix-vector operations, Sn -a±, F ■ (a]")*, and -a~[ . These operations can be performed 
mathematically, but the measurement corresponding to these operations in a physical setting 
is not obvious. Here, we dwell on mapping these matrix-vector operations into their physical 
counterparts, thus making our algorithms physically realizable. 

The first operation, Sn -af, can be realized by measuring the backscattered wave. In an ex- 
perimental setting, the modal coefficient vector of the backscattered wave would be extracted 
from the backscatter intensity measurement by digital holography techniques described in, 
for example [28]. We also assume that it is possible to modulate the amplitude and phase 
of a wavefront, using the methods described in [29]. Thus, the matrix- vector multiplicative 
operation Sn -af corresponds to sending an incident wavefront with modal coefficient vector 
a~l and measuring the modal coefficient vector of the backscattered wavefront. Furthermore, 
we assume that these modal coefficient vectors can be recovered perfectly, and the ampli- 
tude and the phase can be perfectly modulated, so that we might investigate the best-case 
performance of the algorithms. 

The second operation, F ■ (a]~)*, can be realized by time-reversing the wave. Let flipud(-) 
represent the operation of flipping a vector or a matrix argument upside down so that the 
first row becomes the last row and so on, and let * denote complex conjugation. We define 
F = flipud(J), where / is the identity matrix; then the operation F ■ (a^)* represents time- 
reversing the wave corresponding to . This can be explained as follows. The expression for 
time-reversed wave of is 

/ N \*N 

Mp)r= E ^v*""' 2 = E KKnYe^- 



yn=-N / n=-N 

N 



E K{al_ n ye-^-. (14) 



n=-N 

Note that we have used the fact that h*_ n = h n and k_Z n = —k£- From Eq. (14), we 

see that the modal coefficient vector representation of the time-reversed wave of aj~ is 

1 T 

{a^y (a^v-i)* • • • ( a -iv+i)* ( a -7v)* = F ' (&)*■ Furthermore, we emphasize that 
the operation F ■ (aj")* can be physically realized via phase-conjugate mirroring (PCM) [22]. 

The third operation, S^\ ■ a±, can be realized by using reciprocity. In a scattering medium 
that exhibits reciprocity, there are relationships [30-34] between the incident and scattered 
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wavefronts. Consequently, reciprocity requires the reflection matrix Su to satisfy 

S» = F- S* u ■ F (15) 

This means that if a is an input to the system that produces a backscattered wave of b, then 
sending F ■ (a)* will produce backscattered wave of F ■ (b)* in a medium whose reflection 
matrix corresponds to S^. (Fig. 3) 

An important implication of this equation is that the matrix- vector operation ■ aj~ can 
be cast in terms of physically realizable operations. Note that ■ can be expressed as 

S* ■ ar = F ■ S* u ■ F • ar = F ■ (S u ■ (F ■ (&)*))*. 

From the last expression, we see that the operation Sf{ ■ % can be physically realized in a 
sequence of two steps: 

1. Time-reverse the wavefront whose modal coefficient vector is a~[ , and send it to the 
scattering system. 

2. Time-reverse the resulting backscattered wavefront. 

We call this sequence of operations as double phase conjugation, and we shall leverage it 
extensively in what follows. 

5. Iterative, physically realizable algorithms for transmission maximization 

We now develop iterative, physically realizable algorithms for transmission maximization that 
converge to a opt in Eq. (9), by utilizing significantly fewer backscatter field measurements 
than the O(M) measurements it would take to first estimate Su and then compute its SVD 
to determine v_ M . 

5. A. Steepest descent method 

The backscatter minimization problem involves optimization with respect to the objective 
function ||Sii • af||| that appears on the right hand side of Eq. (8). The objective function's 
negative gradient is used as a search direction to correct the previous input as 

d\\Sn-ai\\t 



^i,(fe+i) — -i,(fc) t 1 



^t(fc) _ 2 /^n " " fi£ 



.(*)' 

-1 -l,(fe) 



where cl[,q represents the modal coefficient vector of the wavefront produced at the k- 
th iteration of the algorithm and fi is a positive stepsize. This yields Algorithm 1 which 
iteratively refines the wavefront a^^^ until the backscattered intensity ||Sii •<2y f ( - fc )||2 drops 
below a preset threshold e. 
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Algorithm 1 Steepest descent algorithm for finding a, 



1: Input: ( ) = Initial random vector with unit norm 

2: Input: fi > = step size 

3: Input: e = Termination condition 

4: k = 

5: while \\S n • a^-Jl! > e do 

6: a+ (fe) = a+ (fc) - 2(iSf[ ■ S u ■ a+ (fc) 

7: ai,(fe+i) = ^i,(fc)/IK,(fc)ll2 
8: k = k + 1 
9: end while 



Armed with the relationship in Eq. (15), step 6 in Algorithm 1 can be expressed as 

oJ( fc ) = fij(fc) - 2 /^ii • Su • = aj( fc ) - ■ ■ F ■ S u ■ a+ (jfc) . (16) 

This allows us to recast each step of Algorithms 1 into the counterparts of the physical 
operations in the second column of Table 1. 

The sequence of steps 1 — 4 in Table 1, which involves double phase conjugation, amplifies 
the highly-backscattering component in the wavefront, analogous to the operations for time- 
reversal focusing [22,35-37]. In step 5, this component is subtracted leading to a refined 
wavefront that will backscatter less. This process is repeated till convergence. A consequence 
of this technique is that the backscatter field intensity will typically decrease monotonically. 
This makes the measurement of the backscatter modal coefficient vector increasingly difficult 
as the iteration progresses. An additional disadvantage of this method is the obvious need 
to carefully set \x to guarantee convergence, < ji < ~iK$ ~ 1. In an experimental setting, 
the step size \i is chosen by a simple line search, i.e., by scanning a set of discretized values 
and selecting the one that results in the smallest backscatter intensity after a fixed number 
of iterations. 

We describe a method next, which maintains high backscatter field intensity throughout 
the process and does not require selection of any other auxiliary parameters to guarantee 
convergence. 

5.B. Conjugate gradient method 

Consider an iterative solution to Eq. (8) where the iterate (before normalization for power) 
is formed as 

^l,(fc+l) = ^t(k) + ^(*+i)^(*)> ( 17 ) 
where fi(k+i) is a stepsize and dn^ is the search direction. In this framework, Algorithm 1 
results from setting fi(k+i) = ^ and d^ = —25^ ■ Su ■ g£r k y 



9 



Vector Operation Physical Operation 



q _|_ -, _|_ Backscatter 

Oj = • a 1(fc) 1 : a 1(fc) ► Oj 

PCM 



1 

2 : = F • (a7/)* 2 : — > 



_i_ Backscatter 

6 : a x = on • o : a-i > a x 



4 : g+ = F • (g7/)* 4 : a+ 

5 : a] 1 " = a|m — 2fia± 5 : = aj^ — 2fiaf 

^ _|_ ~ _|_ / 1 1 ~ _L. 1 1 n ~_|_ Normalization _|_ 

6 : a+ (fe+1) = a+/ ||a+|| 2 6 : a+ ► aj (fc+1) 



Table 1. Steepest descent algorithm for transmission maximization. The first 
column represents vector operations in Algorithm 1. The second column rep- 
resents the physical (or experimental) counterpart. The operation i — > 
F ' {<h)* can be realized via the use of a phase-conjugating mirror (PCM). 
The algorithm terminates when the backscatter intensity falls below a preset 
threshold e. 

The conjugate gradients method (see [26, Chapter 5] for a detailed derivation) results from 
choosing the stepsize 

0(fc+i) = ||£(fc)||!/||Sii ■ d {k) \\l, (18a) 
with the search direction given by 

d( k +i) = L(k+i) + /?(fc+i)d(fe), (18b) 

and 

P(k+i) = ||£(fc+i)|ll/||E(fc)|l2- ( 18c ) 

Here, the residual vector is 



r 



■ S u ■ atr k+1 y (18d) 



(fc+i) — — ^i,(fc+i) 
The iteration terminates when ||r(fc+i)||2 < e, a preset threshold. 

Plugging Eq. (17) into Eq. (18d) and substituting the expressions in Eqs.(18a) - (18c) 
gives us an alternate expression for the residual vector 

£(*+!) = L( k ) - V(k+i)Su ■ S n ■ d {k) , (19a) 

or, equivalently 

1 1 ^* fe) 1 1 ^ 

H(fe+i) = L(k) - Yq 1 — H2^n ' ^ n ' d-(k)- ( 19b ) 

IPll ' «(fc)||2 
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The utility of Eq. (19b) will become apparent shortly. 

To summarize: we described an iterative method for refining the wavefront <Li(k) v ^ a Eq- 
(17). Inspection of the update Eqs. (18a)-(18c) and Eq. (19b) reveals that matrix-vector 
operations Sn ■ d, k \ appears in Eq. (18a) while S X1 ■ S±\ • d, k \ appears in Eq. (19b). This 
means that the vector dr k \ is transmitted and the associated backscatter is measured. Note 
that these measurements are used to iteratively refine the vector , but a|m is never 

actually transmitted until the termination condition | jzirfc+i) 1 1 2 < e is met. This is reflected 
in the physical description of the proposed algorithm in Table 2. Also, note that we start 
with a random unit vector g^^y and set d^ and to —S^ ■ Sn ■ Ow ), since we are using 
conjugate gradient for finding the input that minimizes reflection, i.e., 

_|_ Backscatter _ PCM _|_ Backscatter _ PCM , 

— ^1,(0) * 9li ► Qli ► Qli ► «(o) = £(o)- 



Vector Operation Physical Operation 



1 : 


d_i — Sn ■ d( fc ) 




1 : 


j Backscatter * 

d(k) > «i 


2 : 


d+ = F ■ (a7/)* 




2 : 


,_ PCM ,+ 

d x > dj 


3 : 


d_\ = Sn ' d^ 




3 : 


y_|_ Backscatter 

dj > «i 


4 : 


d = F- (d^y 




4 : 


,_ PCM , 
dry > d 


5 : 


A*(fc+i) - \\L(k)\\ 


l/(dfyd) 


5 : 


^(fc+i) = llH(fc) II l/{df k) -d) 


6 : 


L {k +i) = L( k ) ~ 


H{k+l)d 


6 : 


r {k +i) = L(k) ~ V(k+i)d 


7 : 


P{k+i) = \\r { k + i)h/\\r {k )h 


7 : 


P(k+i) = ll£(*+i)IMl£(fc)l|2 


8 : 


d(k+i) - L(k+i) 


+ /3(fc+i)C?(fc) 


8 : 


d(k+i) = + P(k+i)d(k) 



Table 2. Conjugate gradient algorithm for transmission maximization. The 
first column represents iterates of the conjugate gradients method. The second 
column represents the physical (or experimental) counterpart. The operation 
1 — > F ■ (aj*)* can be realized via the use of a phase-conjugating mirror 
(PCM). The algorithm terminates when the residual vector ||zi(fc+i)l|2 < e, a 
preset threshold at which point the optimal backscatter minimizing wavefront 
is constructed as aw t+1 ] = aw^+A^fc+i)^) followed by a power normalization 

= -t(fc+i)/ll-i,(A;+i)ll 2- 
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A feature of the conjugate gradient method is that the intensity of the backscatter measure- 
ment Sn ■ dr k \ is expected to remain relatively high (for a strongly backscattering medium) 
throughout the process. It is only when the wavefront corresponding to ajf k+1 \ is excited that 
a strong transmission (with minimized backscatter) is obtained - this might be a desirable 
feature for communication or covert sensing applications. Consequently, the algorithm will 
produce high intensity backscatter measurements, thereby facilitating accurate estimation 
of the backscatter modal coefficient vectors that are an important component of the pro- 
posed algorithm. This makes the conjugate gradient method less susceptible to measurement 
noise than the steepest descent method where the backscatter intensity decreases with every 
iteration. 

6. An iterative, physically realizable focusing algorithm 

We first describe a generalized coordinate descent method for amplitude and phase opti- 
mization. Assume we are given a M x Nb matrix B = \b_ 1 ... h.N B ] whose columns are 
orthonormal so that B H ■ B = I^ B . Thus Nb denotes the number of (orthonormal) bases 
vectors. 

The key idea here is to expand af on the right hand side of Eq. (11) in terms of the bases 
vectors given by the columns of B as 

N B 

ai = Y,Pi^ l % (20) 
i=i 

where pi > and 4>i G [— tt, tt] are the unknown amplitudes and phases, respectively. 

The optimal amplitudes can be estimated by transmitting af = 6 ; for every 1 = 1,... Nb, 
measuring the corresponding intensity X\ at the target, and setting pi = a/Z^. This can be 
accomplished with O(Nb) measurements. 

The phases can be estimated by first setting (j>\, ■ ■ ■ <j>N B randomly and then for I = 
1,...,Nb, sequentially finding the phase that optimizes measured intensity. This can be 
done via a simple line search, i.e., by scanning the measured intensity over a fixed set of 
discretized values of the phase or by using more sophisticated algorithms such as golden 
section search algorithm with parabolic interpolation [38, Section 10.2]. This too requires 
O(Nb) measurements. 

Setting Nb = M and B = I yields the coordinate descent approach used by Vellekoop and 
Mosk [6,7,25]. This corresponds to exciting one plane wave mode at a time and inferring the 
optimal phase and amplitude one mode at time. Such an algorithm requires O(M) iterations 
to yield the optimal focussing wavefront. Setting B to the 2D Walsh function basis matrix 
yields the method proposed by Stockbridge et al. in [18]. 

An important insight from Eq. (13) is that if we were to express the optimal focusing 
wavefront as a superposition of eigen-wavefronts of S21, then typically only K <C M of the 
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combining coefficients will be large. Thus only K of the pi coefficients in Eq. (20) will be 
significant if we set B to be the right singular vectors of S21. In the lossless setting, the K 
eigen-wavefronts of S21 that have transmission coefficients close to one correspond precisely 
to the K eigen-wavefronts associated with Su that have reflection coefficients close to zero. 
Hence, we can set B to be the right singular vectors of Su and expect only K of the pi 
coefficients in Eq. (20) to be significant as well. Thus, we need to measure the K singular 
vectors of Su associated with its K smallest singular values. 

The Lanczos algorithm is an iterative algorithm for accomplishing just that [26,27]. The 
key idea is to create a tridiagonal matrix H whose eigenvalues and eigenvectors (referred to 
as the Ritz values and vectors) are approximations of the eigenvalues and eigenvectors of 
Six ■ S\\. The algorithm is summarized in the first column of Table 3; its physical counterpart 
is described in the second column. The matrix B in Eq. (20) is obtained as 



where Q 



and U 



«(1) • • • U(N B ) 



B = Q-U, (21) 

are the Nb vectors produced by the algorithm (see Table 3) 
are the Nb eigenvectors of H associated with the Nb smallest 



eigenvalues. 

The convergence theory [27] of the Lanczos algorithms predicts that the eigenvector es- 
timates will rapidly converge to the K eigenvectors of S^ ■ Su associated with the eigen- 
wavefronts of Su with the smallest reflection coefficients; hence, setting Nb = O(K) will 
suffice. An estimate of K can be formed from the eigenvalues of H by counting how many 
of the converged eigenvalues of H are below a preset threshold e. 

Estimating these K right singular vectors will require O(K) measurements and when K <C 
M, we shall obtain a near-optimal focusing wavefront using significantly fewer measurements 
than the O(M) measurements required by the coordinate descent when B = I. We shall 
corroborate this convergence behavior using numerical simulations next. 

7. Numerical simulations and validation of the existence of highly transmitting 
eigen-wavefronts 

To validate the proposed algorithms, we compute the scattering matrices in Eq. (2) via 
a spectrally accurate integral equation solver that characterizes fields scattered from each 
cylinder in terms of their traces expanded in series of azimuthal harmonics, and that uses a 
fast, Shank's transform-accelerated summation technique to efficiently evaluate all periodic 
Green's functions. In the numerical experiments below, care was taken to ensure 11-th digit 
accuracy in the entries of the computed scattering matrices. 

Fig. 4 shows the empirical transmission coefficient distribution, i.e., the singular value 
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Table 3. The Lanzcos algorithm and its physical counterpart which computes 
a tridiagonal matrix H whose eigenvalues and eigenvectors are closely related 
to the eigenvalues and eigenvectors of ■ Sn. Note that we initialize the 
algorithm by setting k = 1, q.^ to a random unit norm vector, and S(o) = 0. 

squared of the S 2l matrix of a slab with D = 197A, L — 1.2 x 10 4 A, r = 0.1 1A, N c = 14, 000 
(Dielectric), = 1.3, M = 395 and , I = 6.7A, where I is the mean of the minimum- 
inter-scatterer-distances. The computation validates the bimodal shape of the theoretical 
distribution in Fig. 2. 

Next, we consider scattering system with D = 14A, L = 5.4A,r = 0.11A,iV c = 50 (PEC), 
M = 27, and I = 0.8A. Here r norma | = 0.49 while r opt = 0.9995 so that wavefront optimization 
produces a two-fold increase in transmissited power. Fig. 5a and Fig. 5b show the wavefield 
produced by the optimal wavefront and a normally incident wavefront, respectively. Fig. 6 
shows the modal coefficients of the optimal wavefront corresponding to Fig. 5b. 

Fig. 7 displays the rate of convergence of the algorithm's developed for a setting with 
D = 197A,L = 3.4 x 10 5 A,r = 0.11A,JV C = 430,000 (Dielectric), n d = 1.3, M = 395 and, 
I = 6.69A; this slab has a comparable (slightly lower) packing density than that in Fig. 5a. 

A normally incident wavefront results in a transmission of r norma | = 0.038. The optimal 
wavefront yields r opt = 0.9973 corresponding to a 26-fold increase in transmission. Algorithms 
1 and 2 produce wavefronts that converge to the near optimum in about 5 — 10 iterations, 
as shown in Fig. 7. 
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Fig. 8 plots the transmitted power after the 10-th iteration of Algorithm 1 for different 
choices of [i. Fig. 8 reveals that there is broad range of \i for which the algorithm converges 
in a handful of iterations. We have found that setting \i « 0.5 yields fast convergence. 

The conjugate gradient method (Algorithm 2) converges slightly faster than the steepest 
descent method (Algorithm 1) in the setting where we chose the optimal fi = 0.5037 for 
Algorithm 1 by a line search; i.e., we ran Algorithm 1 over a fixed set of discretized values 
of \x between and 1, and chose the optimal \x that gives the fastest convergence result. 
In an experimental setting, the line search for finding the optimal fi for the steepest de- 
scent algorithm will require additional measurements. Thus, Algorithm 2 will require fewer 
measurements than Algorithm 1 with the additional advantage of not requiring any auxiliary 
parameters to be set. 

Fig. 9 considers the same experimental setting as in Fig. 7 with a target at (D/2, 5.4A) and 
plots the focus achieved at the target by exciting a focusing wavefront as in (12). The modal 
coefficients are plotted in Fig. 10a. Fig. 10b shows the sparsity of the modal coefficients of the 
optimal focusing wavefront when expressed in terms of the basis given by the right singular 
vectors of the Sxx matrix (or equivalently, the eigenvectors of S^ ■ Sxx)- 

Fig. 11 plots the intensity achieved when using Nb bases vectors for the algorithms de- 
scribed in Section 6 in the same experimental setting as in Fig. 9. The new algorithm which 
utilizes the eigenvectors of S^-Sxx associated with its smallest eigenvalues reaches the optimal 
intensity with significantly fewer iterations than the coordinate descent algorithm. Furthe- 
more, the number of iterations taken will not depend on the number of modes controlled 
unlike coordinate descent methods that use fixed bases vectors. 

Finally, we consider the setting where the scatterers are absorptive. Here, backscatter min- 
imization as a general principle for increasing transmission is clearly sub-optimal since an in- 
put with significant absorption can also minimize backscatter. We defined gain as r opt /r norma |. 
Here we have D = 197A,L = 3.4 x 10 5 A,r = 0.11A,iV c = 4.3 x 10 5 (Absorbing Dielectric), 
rid = 1-3 — jn, M = 395, and I = 6.69A. In Fig. 12, we compare the gain obtained by using 
the backscatter minimizing wavefront to the gain obtained by the optimal wavefront (that 
utilizes information from the S21 matrix) for various n, as the thickness of the scattering 
system increases. We obtain an increase in transmission and the methods described again 
produce dramatic gains whenever the scatterers are weakly absorptive. 

8. Conclusions 

We have numerically verified the existence of eigen-wavefronts with transmission coefficients 
approaching one in highly scattering systems and developed physically realizable algorithms 
for finding these highly transmitting eigen-wavefronts using backscatter analysis. We also 
developed a physically realizable algorithm for forming a focused input using the highly 
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transmitting eigen-wavefronts identified by the previous algorithm. Via numerical simula- 
tions it was shown that the algorithms converged to a near-optimal wavefront in just a 
few iterations. The proposed algorithms are quite general and may be applied to scattering 
problems beyond the 2-D setup described in the simulations. We are currently investigating 
extensions to imaging and sensing applications. 
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Fig. 1. Geometry of the scattering system considered. 
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Fig. 2. Theoretical distribution in (7) for L/l = 3. 
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Fig. 3. The relationship between wavefronts in a medium that exhibits reci- 
procity. Reciprocity tells us that ■ a is obtained by time-reversing the wave 
before and after sending a into the medium, and we call this sequence of op- 
erations double phase conjugation. 
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Fig. 4. Empirical transmission coefficients distribution from a scattering system 
with D = 197A,L = 1.2 x 10 4 A,r = 0.11A,iV c = 14,000 (Dielectric) ,n d = 
1.3, M = 395,7 = 6.7A, where I is the mean of the minimum-inter-scatterer- 
distances. 



20 



15 



10 



X 







□ 

9- 


>■ 




r 




.' 


... , • - 

• 

• A* 


o 


* % 
















. % 1 


o »*° • 


■ ■ 




, • • •> 


ill ] 
4 1 


> ° . • 

f o ° 
ft o 

OA. 


c 








□ 














-5 





5 10 



(a) Wavefield produced by a normally incident wavefront. 
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(b) Wavefield produced by the optimal wavefront. 

Fig. 5. Wavefield plot of the incident-plus-backscatter wave corresponding to 
(a) normally incident and the (b) optimal wavefront, which were sent to a 
scattering system with D = 14A, L = 5.4A,r = 0.11A,7V C = 50 PEC, M = 
27,7 = 0.8A. The normally incident wavefront has r norma | = 0.49 while the 
optimal wavefront yields r opt = 0.9995. 
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Fig. 6. The modal coefficients of the optimal wavefront corresponding to Fig 
5b are shown. 
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Algorithm 2 (Conjugate Gradient) 
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Fig. 7. The transmitted power versus the number of iterations is shown 
for steepest descent algorithm with fi = 0.5037 and for conjugate gradi- 
ent in the setting with D = 197A,L = 3.4 x 10 5 A,r = 0.11A,iV c = 
430, 000 dielectric cylinders with nrf = 1.3, M = 395,7 = 6.69A. The conju- 
gate gradient algorithm converged to the optimal transmitted power slightly 
faster than the steepest descent algorithm. However, since the steepest de- 
scent algorithm requires a line search for setting the optimal step size /x, it 
requires more measurements than the conjugate gradient method which does 
not require any parameters to be set. 
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Fig. 9. Intensity plot around the target at (D/2, 5.4A) outside the scattering 
system defined in Fig. 7. The optimal focusing wavefront forms a sharp focus 
of 1A around the target. The unoptimized wavefront solution corresponds to 
an incident wavefront that would have produced a focus at the target if there 
were no intervening scattering medium. 
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Fig. 10. Here, we depict the magnitude of the coefficients of the optimal focus- 
ing wavefront, corresponding to the situation in Fig. 9, in terms of two choices 
of bases vectors. In (a) we decompose the optimal focusing wavefront with re- 
spect to the bases vectors corresponding to plane waves; in (b) decompose the 
optimal focusing wavefront with respect to the bases vectors associated with 
the eigen-wavefronts of the Sn matrix. A particular important observation is 
that the eigen-wavefront decomposition yields a sparse representation of the 
optimal focusing wavefront. 




Fig. 11. Intensity at target as a function of the number of bases vectors for 
the new algorithm (which uses the bases vectors estimated using algorithm de- 
scribed in Table 3) versus the standard coordinate descent method which uses 
the plane wave associated bases vectors (see Section 6) for the same setting as 
in Fig. 9. The sparsity of the optimal wavefront's modal coefficient vector when 
expressed using the bases of the eigen-wavefronts (shown in Fig. 10b) leads to 
the rapid convergence observed. The optimal wavefront was constructed as 
described in Section 3.B using time- reversal. 
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Fig. 12. Gain (=:r opt /r norma |) versus the thickness L/X in a setting with 
D = 197A, r = 0.11A, N c = 430, 000 Absorbing Dielectric, n d = 1.3 - jn, M = 
395,1 = 6.69A, for different values of k. The solid line represents the maxi- 
mum possible gain and the dashed line represents the gain obtained by using 
backscatter minimizing algorithm discussed in Section 5. 
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